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Abstract 

Many simulations of stochastic processes require colored noises: I describe here an exact numer- 
ical method to simulate power-law noises: the method can be extended to more general colored 
noises, and is exact for all time steps, even when they are unevenly spaced (as may often happen 
for astronomical data, see e.g. N. R. Lomb, Astrophys. Space Sci. 39, 447 (1976)). The algorithm 
has a well-behaved computational complexity, it produces a nearly perfect Gaussian noise, and its 
computational efficiency depends on the required degree of noise Gaussianity. 
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In recent years colored noise sources have been considered in many disparate applications, 
that range from stochastic resonance [lj , to biophysics and beam dynamics in particle 
accelerators 0, 0]. The analytical approach to some of these processes is often difficult, 
and sometimes impossible, and numerical experiments are commonly used to support the 
analytical conclusions, or as an aid to discover new results. For this reason, algorithms 
that produce colored noise have acquired an ever increasing importance. This widespread 
interest spans different scientific communities, and the existing algorithms reflect the vari- 
ety of approaches to the understanding of stochastic processes in different contexts. There 
are physics-inspired algorithms that rely mostly on equations of the Langevin type, FFT- 
based and autocorrelation function methods that use the spectral or correlation properties 
of colored noise, and time-series methods that produce colored noise from different filtering 
approaches. The review paper by Kasdin 6] provides a long list of references until 1995, 
centered mostly on linear processes and FFT methods. More recently, Greenhall wrote a 
review paper on FFT-based methods Q| , and reference js| is another very clear paper on the 
same topic. I describe here an exact numerical simulation of power-law noises that can be 
extended to more general colored noises, and which is based on the classical argument pro- 
posed longago by Bernamont to model l/f a noise as a superposition of Ornstein-Uhlenbeck 
processes |9|. The synthesis of colored noise from a point process is clearly not new, because 
this kind of modeling dates as far back as 1909, to the work of Campbell (see also the 
famous paper by Rice 3])! more recently Teich, Lowen and collaborators have carried out 
extensive studies on point processes with long-tail pulse response functions |l^, Q|, and 
others have studied the synthesis of power-law spectra from nonlinear processes (see, e.g., 



c 16 for a model based on a multiplicative point process). The simulation methods described 
in (| [?| assume evenly distributed sampling steps, and the extension to uneven sampling is 
not trivial: however noneven s ampling has many important applications (see, e.g. the classic 
papers by Lo m b and Seargle HQ on period ana,y S is for nregnlarly sampled astrononr- 
ical data, and two more recent references |19j,|20j), and Gillespie discussed a method valid 
for the Ornstein-Uhlenbeck process and based on the Langevin equation in 1996 The 
algorithm proposed here is very general (it is not limited to the OU process), it is very easy 
to implement, it is valid for all time steps, it has a well-behaved computational complexity, 
and produces a nearly perfect Gaussian noise, and its computational efficiency depends on 
the required degree of noise Gaussianity. 
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Here we take a signal x(t) that originates from the linear superposition of many random 
pulses, i.e., pulses that are random in time and can be described by a memoryless process 
with a Poisson distribution, have random amplitude A drawn from a distribution with finite 
variance and probability density g A (A), and such that their pulse response function is 

, | exp(-At) if t > 

h(t,X) = { ~ (1) 

j if t < 

with a decay rate which is drawn from a distribution with probability density g\(X), so that 

x(t) = J2 A kh(t-t k ,X k ) (2) 

k 

where t k is the time at which the /c-th pulse occurs, A k is its amplitude and X k is the decay 
rate of its pulse response function. 

If n is the pulse rate, then on average there are n [gA{A)dA] [g\(X)dX] dt pulses in the 
time interval (t', f + dt) and in the amplitude-A range dAdX; the number of pulses follows 
a Poisson distribution and therefore the variance of the number of detected pulses is also 
equal to n \g^{A)dA] [g\(\)d\] dt. This means that the mean and the variance of the output 
signal at time t are given by the integrals 

P^max pA ma ,x pt 

(x) = g\(X)d\ / g A (A)dA dt'n[Ah(t -t',X)\ (3) 

and 

"j4.m.ax pt 

n2 



((Axf) = gx{X)d\ g A (A)dA dt'n[Ah(t - 1' , A)f (4) 

" ^min J A m i n J OO 

If we assume that the amplitude A is fixed, we take the pulse response function (Q), and 
rearrange the time integration, the integrals (jUJ) and (jlj) simplify to 



^max pOO P^ r , 

gx 

A "~"\A 



(x)=nAl g\(X)dX I dt [h(t, A)] = nA I 9^dX = nA{\) (5) 



and 



r^max POO pXmax „ I \ \ „ A2 / 1 \ 

((Ax) 2 ) = nA 2 / g x (X)dX / dt [h(t, A)] 2 = nA 2 / ^dX — I -J (6) 

Now let H(u, X) be the Fourier transform of h(t, A), then from the causality constraint on 
h(t, X) and Parseval's theorem we find that the variance (jUJ) can be trasformed into 

„/l2 rXmax POO a2 POO rXmax 

((Ax) 2 ) = —I g x (X)dX I dw\H(u,X)f = ^ / dw \ g\(X)dX \H(lu, X)\ 2 

(7) 
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The right-hand expression in equation (JJJ) shows that the spectral density is 

SM = — / g x (X)dX\H(uj,X)\ 2 (8) 

and since \H(u, A)| 2 = (u 2 + X 2 )~ 1 for the exponential pulse response function ([U), we obtain 
eventually 

We consider now three special, important cases: if there is just a single decay rate A, the 
spectral density has the usual Lorentzian shape 

n A 2 1 

s ^ = (10) 

which, for u ^> A has a l// 2 behavior, so that we can approximate a l// 2 spectrum in an 
actual process by choosing a A smaller than the lowest observed frequency. With a careful 
choice of the distribution g\(X) we can synthesize many different spectra, but there are two 
special choices for g\(X): we can take a uniform distribution or a range-limited power-law 
distribution. If we assume a uniform distribution of decay rates, between X min and X max , 
i.e. 

9x(X) = - x ^— - (11) 



then the average (1/A) that determines the mean level (jSJ), and the variance © is 

1 \ ln( A max / A m j n ) 



A / X max A m j n 



(12) 



and using equation El the spectral density is easily shown to be 

tlA? 1 / A mm A 



ol \ j. "max , "mm i / 1Q i 

o{u!) = r— I arctan arctan I (13) 

^{X max — Xmin) LU \ LO LO J 

and in the range X min <C io <C X max this spectral density has a 1/ f behavior. Similarly, if 
we take a range-limited power-law distribution 



^)=[ 'o-j A-' (11) 

Amax *min 



then the average (1/A) is 

1 \ /I — /?\ A - ' 3 — A - ' 3 

\ / \ ^max min 



^ I \ ft J Xmax A • 



(15) 
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and the spectral density is 
S(u>) 



\l-/3 p I * l± -i 
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(16) 



""" v 2 ' ' 2 ' cu 2 
which has a 1// 1+/3 behavior in the range \ m in "Cw-C A maa; . 

Now we follow the lead provided by these considerations, and we take, e.g., the case where 
there is just a single decay rate A, so that the spectral density has the Lorentzian shape (JTUJ). 
Since the probability density of the time intervals Atk between Poisson events is well-known 
to be dP(At) = nexp(—nAt)dAt, we can generate a sequence of At's from an exponential 
distribution, and we can thus generate the sequence {tk} (with tk+i = tk + Atk) required to 
evaluate a realization of x(t) as in equation (J2J): figure shows an example where the single 
decays are clearly visible. Figure [T] also shows that, although the process has the desired 
spectral density, it is quite obviously non-Gaussian and therefore this generation method 
seems to be of marginal utility, as most of the actual physical processes are Gaussian and 
Gaussianity is usually a required property of a good noise generator (see, e.g., the recent 



paper 



111 ] that describes a hardware-based Gaussian white noise simulator and contains a 
list of relevant references: notice also that Gaussianity is sometimes a weakness rather than 
a strength, see, e.g. [l2|). The Gaussianity in shot noise processes has been studied at 
length since the paper by Rice and here we strictly limit the discussion to the special 
processes considered in this paper. The single exponential spikes in figure Q stand out more 
clearly when the average rate n of the Poisson process is smaller than the decay rate A; by 
contrast, when n ^> A, at any time there are many pulses of comparable size and the sum 
has a nearly Gaussian behavior. We can gain further insight in this generation method by 
using the mean moment generating function (mmgf) for a Poisson process with average rate 
a: 

(exp [it{k - (k))}) = J2^ k ~ ( k )) m )^T = ex P [( e " " l ) a ~ lta \ ( 17 ) 

m,=0 

%^ Ij^ %^ 

= 1 + —at 2 + —at 3 + —(3a 2 + a)t 4 
2! 3! 4P ; 

+ ^y(10a 2 + a)t 5 + ^(15a 3 + 5a 2 + a)t 6 + 0(t 7 ) (18) 

Now we use the mmgf to compute the higher-order moments: as already discussed in the 
derivation of equations (jUJ) and (jlj), the process x(t) is the sum of Poisson variates with 



different amplitudes; on the other hand the mmgf of the weighted sum ak + j3j of two 
independent Poisson variates k and j, both with rate a, is 



(exp {it [(ak + (3j) - (ak + /3j)]}) = (exp [ita(k - (£;>]> (exp {it(3(j - (j)}) (19) 

= 1 + ~a [(at) 2 + ((3tf] 

+|j [(«t) 3 + ((3tf] + 0(t 4 ) (20) 

Using the mmgf's given above we could proceed as in standard texts on probability theory, 
and show that for large a the process approaches an exact Gaussian distribution (the usual 
proof of the Central Limit Theorem), but the purpose here is giving a quantitative estimate of 
the deviation from Gaussianity: from the expansion (J2(J|) . we see that, just like the variance, 
the third moment about the mean of the weighted sum of independent Poisson variates is 
the weighted sum of the third moments of the individual variates (this is not true for the 
fourth and the higher moments), and therefore we can write a simple expression for the third 
moment about the mean 

((Ax) 3 )=nA 3 9x(X)d\ dt[h(t,\)f (21) 

and we can use this expression to compute the skewness of the frequency distribution 

skewness = « A *> 3 > = fcxw*r*w.*)f (22) 
{{Ax)2)3/2 {it: r * w*. A )i 2 f 2 

With the pulse response function ((TJ) the integrals are easily evaluated, so that 

skewness = ^ Am '" — — ^ = = (23) 

fo\[* ma *[g x {\)/(2\)d\\ 3 V n (VA) 
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From equations ()22|) and (|2*3*j) . we see that the skewness is small when n(l/A) is large (as it 
should be for a Gaussian distribution) and that the actual amount of skewness depends on 
the adimensional product n(l /A), as expected; as a rule of thumb one might take n(l /A) > 10 
for good Gaussianity. 

The previous considerations apply to the noise process x(t) without any reference to 
sampling, however the simulation of noisy physical systems usually implies evaluating the 
noise process at evenly spaced sampling times, so we take now a sequence of sampling times 
{sj} with average sampling interval (As). At each sampling time only the recent pulses 



actually contribute, while the older pulses quickly fade away, for instance the average total 
contribution of pulses that are older than A 'decay / A is 

(5x) = nA 9x(\)d\ / dt'h(t-t',X) 

p^max POO / 1 \ 

= nA 9x(\)d\ / h{t")dt" = nAl-\ e ~ Nd ^ (24) 

which is just a small fraction of the mean value: (5x) / (x) = e~ Ndeca y, and for this reason as we 
proceed forward in time, we can just forget the older transitions. In an actual implementation 
we fix A decay, but because of random event clustering we cannot know a priori how many 
transitions times must actually be kept in memory: for this reason the Poisson-distributed 
transition times should not be stored in an array, but in a linked list [22]; the linked list 
must also store the decay rates that correspond to each transition event. At each sampling 
step the list is updated first by generating as many transition times (and the associated 
decay rates, which are drawn from a given decay rate distribution) as needed to reach (and 
possibly surpass) the actual sampling time Sj, and then by discarding those events with an 
occurrence time t k such that Sj — t k > N decay /\ (see figure |2J). The mean list length is just 
nNdecayi^/ A), and the processing time is proportional to the number of list elements. At 
startup the list is empty, and the first A r decay / ^min(As) samples must be used for initialization 
and afterwards discarded as the algorithm fills the list up to the average level, and thus, 
for a desired number of samples N s , we must generate a total of N s + A decay /(A min (As)) 
samples. The time-complexity of the algorithm is thus proportional to the sum of the total 
number of generated transitions plus the total number of operations used for the list scans, 
i.e., 

complexity = (n s + ^^yj—^j (c^As) + C 2 nN decay ^~ 

N d ecay 1 \ I r~i ^ decay I 1 



= ^l ff " + WuJl a+G MWJ (25) 

The algorithm described above is easily implemented; figures 01 to H show the results 
obtained in a simulation of N s = 2 18 = 262144 transitions, with a single decay rate A = 0.001, 
and a Poisson transition rate n = 1 (here and in all the following discussions the system 
of units is arbitrary); moreover A decay = 20, so that the average relative error due to the 
past transitions that have been discarded is (5x) / (x) = exp(— A decay ) rs 2 • 10~ 9 . With 
these parameters we expect an average list length nN d ecayl A = 20000, and a corresponding 



filling time Ndecay/ (AAs) = 20000: figure 01 shows the list length, which behaves exactly 
as expected. Figures 0] and |S] show the normalized signal amplitude (x(t) — (x))/a, where 
a is the standard deviation of the amplitude, i.e., the square root of the variance (jUJ): at 
the beginning the linked list which contains the process memory is empty, and the signal 
is very far off the predicted average, but as the list fills up to level, the signal quickly 
reaches the predicted average. Figure is the histogram of the normalized signal amplitude 
obtained from 262144 samples, after the list fill-up; the continuous curve superimposed on the 
histogram is a Gaussian with the mean and standard deviation estimated from the samples, 
and we see that there is no visible skewness, because in this simulation run X/n = 0.001, 
which corresponds to a very low skewness (|23|) . but there are multiple peaks, which are 
due to the nonstationarity of a true 1/f 2 process (which is well approximated here), and 
which require an extremely long observation time to establish the Gaussianity of the process 
I23I ] . Finally figures [7| and |H1 show the Discrete Fourier Transform (DFT) spectrum of the 
normalized signal amplitude, which reproduces quite well the expected shape (tTt)|) . 

The next set of figures shows the results of a simulation with a decay rates that are 
uniformly distributed between X m i n = 0.0001 and X max = 1. Just as before, the simulation 
contains N s = 2 18 = 262144 samples, with a transition rate n — 10, a sampling interval 
As = 1 and Ndecay = 20: from these parameters we obtain the average (1/A) ~ 9.211, so 
that the fill-up time is N decay / X min = 200000, the fill-up length is n(N decay (l/ X)) ~ 1842., 
and the number of samples required for the initial fill-up is (Ndecay/ X m i n )/ As = 200000. 
Figure shows the linked list length, which in this case does not reach the average filling 
level with a linear growth law, but with a smoothed curve. Figure El shows the initial part of 
the simulated signal, and figure ^2 shows the histogram of the normalized signal amplitude: 
in contrast to the histogram in figure El now the amplitude distribution is slightly skewed, 
because in this simulation run l/(n(l/X)) ~ 0.109, much higher than the calculated skewness 
for figure H3 Finally the averaged DFT spectrum is shown in figure IT2l the spectrum mimics 
quite well the behavior of a true 1/f spectrum over more than three frequency decades. 

The El shows the average spectrum obtained in a long simulation run with a different 
power-law noise: N s = 2 20 = 1048576 samples have been generated with A = 1 and a range- 
limited power-law distribution of decay rates ()14|) with f3 = 0.2, in the range X m i n = 0.0001, 
^max — 1- Here too the spectral resolution Auj « 0.00038 is larger than the minimum decay 
rate A m j n , and the noise samples reproduce the behavior of a true 1/f 1 ' 2 spectrum (dashed 
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line), over more than three frequency decades. 

This generator can be used to test a standard hypothesis that is commonly used with 
FFT-based colored noise generators, in analogy to the well-known behavior of white noise, 
namely that the standard deviation of the real and imaginary parts of the discrete Fourier 
components Fh of a colored noise process is proportional to the square root of the noise 
spectrum Sk (of. Figure ITU shows the ratio var[3?(Ffc)]/5'fc for a simulated 1/ f noise: the 
average ratio is constant and thus the simulation does not disprove the standard assumption, 
at least for this particular noise process. 

In all the examples described above the sampling interval As is fixed, but the method is 
in no way limited to constant sampling intervals. And indeed this is probably the greatest 
strength of this noise generator, its ability to work also with uneven sampling intervals: this 
is not true for the other common generators 

To conclude, in this paper I have described a generator of colored noise that is exact, is 
not limited to evenly distributed samples, has a well-behaved complexity 0(N S ) (in contrast 
to many other generators that have a 0(N s logN s ) complexity), and is not troubled by 
hidden periodicity issues, like the FFT-based generators. 
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time (arb. units) 



FIG. 1: A realization of the random process x(t) (eq. |2j) with A = 1, n = 1, and with fixed decay 
rate A = 0.5 (all quantities are given in arbitrary units; A is given in inverse time units). The single 
exponential decays are clearly visible, and the random process is obviously non-Gaussian. 
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FIG. 2: Structure and update dynamics of the linked list that holds the Poisson-distributed 
transition events: a. structure of the list at the j-th sampling time Sj] each node contains a 
variable that points to the next node (the end of the list is marked by a null pointer), and stores 
the time tk of the transition and the decay rate of the associated pulse response function. The 
list contains only nodes such that Sj — < Nd ecay /\k. The response of the system is computed 
from the sum exp[— Xk(sj — tk)], where the index k ranges over all the list elements such that 
Sj — tk > (the list head is usually excluded), b. if the (j + l)-th sampling time is greater than 
the time stored in the list head (as is usually the case), the program generates as many transition 
times as needed to reach (and possibly overcome) the (j + l)-th sampling time (light-gray boxes 
in the figure, primed quantities), and next it scans the list to discard all the nodes such that 
Sj+i — tk > Ndecay/^k (dark-gray boxes). At this point the program computes the new response 
and steps to the next sampling step. 
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FIG. 3: Length of the linked list in a simulation with A = 1 and a single decay rate A = 0.001: 
the linked list is initially empty, at it fills up at a constant rate. In this case n = 1, N decay = 20 
and As = 1 and therefore the fill-up time is {Na< eca y/\) = 20000 the fill-up length is niN^ecayl A) = 
20000, and the number of samples required for the initial fill-up is (A^ em j / /A)/As = 20000. After 
the initial fill-up the length of the linked list fluctuates about the average filling level. 
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FIG. 4: Plot of the normalized signal amplitude (x(t) — {x))/a {a is the standard deviation of 
the amplitude, i.e., the square root of the variance © ) in the simulation run described in figure 
Eland in the text. At the beginning the linked list which contains the process memory is empty, 
and the signal is very far off the predicted average; as the list fills up to level, the signal quickly 
reaches the predicted average. 
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FIG. 5: Detail of the normalized signal amplitude shown in the figured just after the list has 
filled up to the average level. This signal displays the lar ge c haracteristic upward an downward 



swings that are well-known in the theory of random walks 
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FIG. 6: The histogram shows the amplitude distribution of 262144 samples from the realization of 
the random process x(t) shown in figure 01 after the list fill-up. The continuous curve is a Gaussian 
with the mean and standard deviation estimated from the samples. There is no visible skewness, 
because in this simulation run X/n = 0.001, which corresponds to a very low skewness (|2H|I. but 
there are multiple peaks, which are due to the nonstationarity of a true l// 2 process (which is 
well approximated here), and which require an extremely long observation time to establish the 
Gaussianity of the process 
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FIG. 7: DFT spectrum obtained from 262144 samples from the realization of the normalized 
signal amplitude (x(t) — (x))/a shown in figure |IJ The continuous curve shows the theoretical 
power spectral density (|lUj) . Because of sampling without low-pass filtering there is some aliasing 
and the DFT spectrum shows a slight upward bend at high frequency. Since the sampling interval 
is As = 1, the Nyquist (angular) frequency is just WNyquist = tt (here and in the following spectra 
time is measured in arbitrary units as in the previous figures, and frequency units are defined 
accordingly). The arrow marks the position of the single decay rate in this simulation A = 0.001. 
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FIG. 8: Averaged DFT spectrum obtained from the same 262144 samples as the spectrum in 
figure [7J split in 32 blocks of 8192 samples each. The continuous curve shows the theoretical power 
spectral density (|lf)[). and now it includes also the first-order correction to aliasing. Because of the 
low-frequency correlation between the blocks (that have been obtained from the same simulation 
record), the average spectrum is a bit higher than expected and the theoretical prediction has been 
globally shifted 20% higher to fit the average spectrum; this artifact is absent in the analysis of 
the whole record (the low-frequency plateau of the spectrum in figure fits the theoretical curve 
exactly as expected). As in figure Q the arrow marks the position of the single decay rate in this 
simulation A = 0.001: because of the shorter record length used for DFT analysis, the frequency 
resolution is poorer here, and the spectrum mimics quite well the behavior of a true 1/f 2 spectrum, 
over about three frequency decades. 
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FIG. 9: Length of the linked list in a simulation with A = 1 and a uniform distribution of decay 
rates in the range \ m in = 0.0001, X ma x = 1: the linked list is initially empty, at it fills up with 
a variable rate that depends on the distribution of decay rates. In this case n = 10, N^ecay = 20 
and As = 1 and (1/A) ~ 9.211, and therefore the fill-up time is N^ecayl '^min = 200000 the fill- 
up length is n(Ndeca y {l/X)) ~ 1842., and the number of samples required for the initial fill-up is 
{Ndecay / Xmin) / lS.s = 200000. After the initial fill-up the length of the linked list fluctuates about 
the average filling level. 
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FIG. 10: Detail of the normalized signal amplitude in the simulation of figure^ just after the list 
has filled up to the average level. 
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FIG. 11: The histogram shows the amplitude distribution of 262144 samples from the realization 
of the random process x{t) shown in figure 1101 after the list fill-up. The continuous curve is a 
Gaussian with the mean and standard deviation estimated from the samples. In contrast to the 
histogram in figure EJ now the amplitude distribution appears slightly skewed, because in this 
simulation run l/(n(l/A)) ~ 0.0109, noticeably higher than the corresponding value for figure|Hl 
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FIG. 12: Averaged DFT spectrum obtained from the 262144 samples in the simulation of figure 
IHJ split in 32 blocks of 8192 samples each. The continuous curve shows the theoretical power 
spectral density which includes also the first-order correction to aliasing. The arrows mark 
the positions of the extreme decay rates \ m in = 0.0001 and A max = 1. The spectral resolution 
Alv k, 0.0015 is larger than the minimum decay rate A m j n , and the spectrum mimics quite well the 
behavior of a true 1/f spectrum (dashed line), over more than three frequency decades. 
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FIG. 13: Averaged DFT spectrum obtained from 2 20 = 1048576 samples with A = 1 and a 
range-limited power-law distribution of decay rates X -13 , with (3 = 0.2, in the range \ m in = 0.0001, 
^max = 1) split in 32 blocks of 32768 samples each. The continuous curve shows the theoretical 
power spectral density ()16|) . which includes also the first-order correction to aliasing. The arrows 
mark the positions of the extreme decay rates \ m in = 0.0001 and A max = 1. The spectral resolution 
Alo ~ 0.00038 is larger than the minimum decay rate X m in, and the spectrum mimics quite well 
the behavior of a true 1/f spectrum (dashed line), over more than three frequency decades. 
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FIG. 14: Ratio vav[^t(F k )]/ S k obtained from 2 20 = 1048576 samples with A = 1 and a uniform 
distribution of decay rates, in the range \ m in = 0.0001, A maa; = 1, averaged over the DFT results 
obtained from 128 blocks of 8192 samples each. The average ratio fluctuates about constant level, 
and this is in line with the usual hypothesis that var[5R(F/ c )] oc Sf. (see, e.g., jsj])). 
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